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SUMMARY 

The numerical solution of exterior problems is typically accomplished 
by introducing an artificial, far field boundary and solving the equations on 
a truncated domain. For hyperbolic systems, boundary conditions at this 
boundary are often derived by imposing a principle of no reflection. However, 
waves with spherical symmetry in gas dynamics satisfy equations where 
incoming and outgoing Riemann variables are coupled. This suggests that 
’natural’ reflections may be important. We propose a reflecting boundary 
condition based on an asymptotic solution of the far field equations. We 
obtain nonlinear energy estimates for the truncated problem and present 
numerical experiments to validate our theory. 
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INTRODUCTION 


Interesting and important problems in gas dynamics are often posed in 
exterior domains. Examples include the explosion of gas bubbles in various 
media and flows external to aircraft. An approach to the numerical solution 
of such problems is to restrict the computational domain to a finite region 
through the introduction of an artificial boundary. For large time compu- 
tations interactions between the solution and the artificial boundary can 
strongly influence the results. The focus of this paper is the development of 
an accurate treatment of these conditions. 

A variety of authors have invoked a principle of no reflection. Notable 
among these are Engquist and Majda [2] who studied the general linear case 
and Hedstrom [8] and Thompson [10] who considered nonlinear hyperbolic 
systems. However, as pointed out by Gustafsson and Kreiss [3], conditions 
satisfied by the exact solution may involve reflections. The current study in- 
voles spherical waves which exhibit coupling between incoming and outgoing 
Riemann variables. One expects this coupling to result in natural reflections 
which should be accounted for in an efficient numerical treatment. Indeed, 
Thompson [10] documents the disappointing performance of nonreflecting 
conditions in such cases. 

An alternate approach is to incorporate the asymptotic behavior of the 
solution in the far field. Conditions for linear problems based on far field 
asymptotics have been successfully employed by Bayliss and Turkel [1] and 
Hariharan and Bayliss [7], Our procedure is to develop approximate solu- 
tions to the appropriate weakly nonlinear initial boundary value problem in 
the region exterior to the computational domain. A condition is thus ob- 
tained which includes appropriate reflections at the computational bound- 



ary. Conditions involving incoming waves generated by inhomogeneities in 
the discarded region have also been proposed by Gustafsson [4,5]. 

The particular equations under consideration are the Euler equations for 
spherically symmetric, isentropic fluid flow; 


dp dz _ 2 z 

dt ^ dr r 7 


( 1 . 1 ) 


dz 

di 

Here p is the density, z 
conditions are 


+ + m = - 2 7r (L2) 

is the momentum and f(p) is the pressure. Initial 


p(r, 0) = po(r) and z(r, 0) = z 0 (r), r > r 0 . (1.3) 

We also assume that the computational boundary is located at r = L (L > 
ro) and that proper conditions at ro are specified. Finally we assume that 

Po(r) = poo and z 0 (r) = 0, r > X. (1.4) 

The plan of this paper is as follows: In section 2 we follow the construc- 
tion presented by Whitham [11] to obtain asymptotic solutions in the far 
field and derive the boundary conditions. Nonlinear energy estimates are 
established for the resulting finite domain problem in section 3. Section 4 
contains a discussion of the numerical treatment of the boundary conditions. 
In section 5 numerical experiments are presented for an idealized weak ex- 
plosion problem. Our technique is shown to yield the correct steady state 
for values of L significantly smaller than those required by the nonreflecting 
conditions. In the final section we propose extensions of our conditions for 
the truly three dimensional case. 



DERIVATION OF ASYMPTOTIC BOUNDARY CONDITIONS 


We find it convenient to work with equations involving Riemann vari- 
ables. They are 



R — - + G(p), 

P 


(2.1) 


S = - - G{ P ). 

P 


(2.2) 

Here 

GW = [ ^&dp. 

J P 


(2.3) 

Then equations (1.1) 

and (1.2) take the form 
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+ - ■ 
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(2.4) 
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dt 

■ 

2 sfnfiz 

pr 

(2.5) 


Here we assume that for r > L 


\fnp) > - p - 

That is, the flow is subsonic in the far field. Then R and S are, respectively, 
the outgoing and incoming Riemann variables. Note that the lower order 
terms couple the equations for the Riemann variables. This means that an 
outgoing wave generates an incoming wave. 

To derive the boundary conditions, we consider the initial boundary 


value problem (2.4), (2.5) and (1.4) on the exterior domain r > L with 
boundary condition 


R{L,t) = g(t). 

(2.6) 

Solving this problem yields 


S(L,t) = *•[*(•)]. 

(2.7) 
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That is the incoming variable is a functional of the outgoing variable. 

Equations ( 2 . 6 ) and ( 2 . 7 ) represent an exact boundary condition at r = 
L. However, the explicit form of the functioned T is not known in general. 
Therefore, we construct an asymptotic solution of the exterior problem valid 
for L sufficiently large. (Similar constructions for steady state problems can 
be found in Hagstrom and Keller [6].) Consistent with the known far field 
behavior of solutions of the linearized equations we expand R and S as 
follows: 


R{r,i) 

S{r,t) 


Ro + £iM + + .... 


5o + £iM + + .... 

r r z 


( 2 . 8 ) 

( 2 . 9 ) 


We note that 


Ro = — So — G(poo)- (2.10) 


We further assume 

<7(0 = Ro + ^ (2.11) 

for some function H(t). 

Equations for R\ and Si are given by: 


8R\ f R\ + S\ 
~df + \ 2 r 

dSi ( Ri + Si 

~dT + V 2r 


+ 


y/f'M + P(Poo) { — r -^ ) 

y/f'iPo*) ~ P(Pco) — ^- 5 -) 


dR i 
dr 

dSi 

dr 


0 , ( 2 . 12 ) 

0 . ( 2 . 13 ) 


Here 


PM 


vTM 

4 /'(«) ' 


Following Whitham [11] we have retained 1 corrections to the characteristic 


speed to suppress nonuniformities in the expansion as r approaches infinity. 



Note that the source terms are absent at this order, so Ri and S\ are Rie- 
mann invariants of the approximate equations (2.12) and (2.13). Since the 
characteristics corresponding to 5 i all originate at t = 0, we deduce 


Si(r,t) = 0. 


(2.14) 


That is, we have a simple wave. Now equation (2.12) can be solved for 
Ri using the method of characteristics. The differential equations for the 
characteristics are given by 
dt 


dr 


JFUI) + + P(m} 7 1 ] ’ 1 • (2 15) 


Using the fact that R\ is constant along the characteristics we find that 

(2.16) 


t(r;r) = r + 


Ri(r,t(r-r)) = H(r), 

T _ ~ L ~ g ( T ) M ZJ fl(r) ] 

\/ 7 '(/>°°) 


(2.17) 


where we have introduced 


B(r) = 


H(t)[ I + P( Poo ) ] 


VTM 

We remark that this solution may break down where characteristics inter- 
sect, in which case a shock must be fitted in. In order to compute the first 
nonvanishing correction to the incoming variable we consider the equation 
for 52; 

^ - (/fW + ^ = JfMk. (218) 

To solve equation (2.18) we make a change of variables in which 52 is ex- 
pressed as a function of r and r , This yields the above equation in the 
form 


1 dr dr\ dS 2 

y/Fip^o) dt ,r dr) dr 


D(r,r)^ = R^r), (2.19) 
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which we write as 


- N(T y r)D{r,r)^- = N{t,t)Ri(t). 


dr 


( 2 . 20 ) 


Here 


vf i r ri7~\^ ~ B ^ 1 + + ( L + + V))}] 

n{t,t) = ^ /'(/>«>) 7 n , > 


and 


D(t,t ) = 1 + 


( i i_ .fi B .( r Ll 

V + l 1 '+ B MJ TrmJ 

( P-\)H{r ) 


VT{p™)t 

Again we solve equation (2.20) by the method of characteristics. We inte- 
grate from the first charact eristic, r = 0, where S 2 = 0. We will only need 
to know £2 at the boundary r = L . Thus we obtain 

s 2 {l,t) = ( 2 - 21 ) 

jo 


where 


^ = -ND, t(r ;r,L) = L. 


( 2 . 22 ) 


To put equation (2.21) in a more convenient form, we differentiate with 
respect to r: 

r,L))dt 


9S2i d L T ’ T - = Ri(t)N(t,L) + [ Ri(*)- 


di 


dr 


ds. (2.23) 


Now the integral term is of order due to the presence of . Therefore 
we neglect it along with £ contributions to IV. Using f£(£,r) = 1, S = 
So + and Ri{t) = L(R(L,t) - Rq ) we finally obtain 

OS, 


dt 


r = L = 


y^MjRiWt) - £o) 
21 


(2.24) 



Equation (2.24) is the boundary condition we propose. We also note 
that the relationship 


£ - O(p) = + 0(i), 

P T* 


(2.25) 


can be used to derive a number of asymptotically equivalent conditions. For 


example 

dS z^f'ipoo) 

dt Lp 

(2.26) 

and 

dS ( G(p ) - G(p eoVVFM 

dt L 

(2.27) 


ENERGY ESTIMATES 

We now study the problem on the truncated region [r 0 , £], rewriting 
the field equations (1.1) and (1.2) in a convenient form. Moreover, we take 
ro = 0 for definiteness. Thus the problem under consideration is 

% + = ° ■ <3 - i) 

| + 7> + £ [/W1 = 0 • (3 ' 2) 

Initial conditions are 

/>(r,0) = po{r), z(r, 0) = 2 0 (r) r>0 . (3.3) 

Our boundary condition at r — L has the integrated form 

- = G(p) - G( Poo ) + a f'-ds + p f{G{p) - G(poo))ds. (3.4) 

p Jo p Jo 

Corresponding, respectively, to (2.24), (2.26) and (2.27) we have: 


a = ft — 
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a 


a 


0 = 0, 

0 , 0 = 

Jb 


(3.6) 

(3.7) 


In addition we need to introduce a finiteness condition at r = 0 due to 
the singularity of the equations (3.1) and (3.2). This is accomplished by 
demanding 


z(r } £) — > 0 as r — ► 0 . 


(3.8) 


It is difficult to establish the well-posedness of initial-boundary value 
problems for nonlinear hyperbolic systems. We content ourselves with the 
derivation of bounds on the growth of the total energy of the system. The 
(physical) energy density is defined by: 

lz 2 

E = «- + P<P)> (3.9) 

z p 

where the internal energy e satisfies: 


We also define 



(3.10) 


q = WE. 


(3.11) 


Here the gradient is with respect to the variables p and z. Then 

/ lz* . T'iP) z] T 

q 


+ fry- ?) • < 3 - 12 > 


Taking the inner product of (3.1) and (3.2) with q we obtain 

dE 1 d _ 
dt + r*dr* ~ °’ 


where 


$ = [r 2 z(e(p) + ^)] + r 


2 P 2 


(3.13) 
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Now integrating (3.13) over [0,Z] with the weight r 2 we obtain 


*'(0 + ®lo = o, 


(3.14) 


where 


fL 

$(f) = / r 2 E dr, 

Jo 

which is the total energy of the system. It is useful to rewrite $ : 


$ = r — 


Pe + / + 5 - 

2p 


Clearly, (3.8) implies that *£ = 0 at r = 0 leaving us with 
*'(<) = -nL,z(L,t),p(L,t)). 
We now state the main theorem of this section. 


(3.15) 


(3.16) 


(3.17) 


Theorem 1 There exists a bounded function F ( £ ) , depending only on f and 
poo, so that the total energy of any generalized solution of (3.1)-(3.4), (3.8) 
satisfies: 

*(t) < *(0) + T(t). (3.18) 


Proof: 

Using Fubini’s theorem it may easily be verified that (3.4) solved for - p yields: 
= G{p{L,t))-G{p„) + (a + IS)£ e « t -'\G{p{L,3))-G{p 00 ))ds. 

Since the bracketed quantity in (3.16) is positive, the right hand side of 
(3.17) can be positive only if - p is negative. As p is nonnegative and G is a 
nondecreasing function, a positive contribution to the total energy implies: 


G{p) < <?(/>«,)(! + A(0), 


P < G~ l (G(poo)(l + A(<))) , 



and 


z 

P 


< G( Poo )(l + A(<)), 


where 

^(e at -l), 

/?£, a = 0 . 

Substituting these inequalities into (3.17) and integrating from 0 to t imme- 
diately yields the desired result. If we specialize to the case of a poly tropic 
gas, f(p) = kp y , we may further conclude that T grows algebraically if (2.27) 
is employed and exponentially in the other cases. 



NUMERICAL PROCEDURE 


In this section we briefly discuss a particular numerical implementation 
of the boundary conditions we have developed. We note that many different, 
stable implementations are possible. Those we present here are used in the 
numerical experiments which follow. 

Introducing a uniform mesh: 


Ti = (i — 1) • Ar, i = 1,...,JV + 1, 
we denote our approximate solution vector by: 

u! = 


p(r, 

z(r, 


«*) ) ’ 


(4.1) 


(4.2) 


and also introduce notation for the approximate fluxes and sources: 


F(U) = 

U +*/('))■ 

(4.3) 


( \ 


H(U,r) 

II 

(4.4) 


U 



We employ the two step Lax-Wendroff method (Sod [9]) in the interior: 

= W + U Ui) - %-AnUh a) - F(Uf)) 

(4.5) 

U* +At = Uf - ^ T {F{U^f ) - F{U^_f )) 

T 2 2 

(4.6) 

-A tH + U*+f),r^,i = 2,...,N. 

This is second order in space and time for smooth solutions. The boundary 
conditions are used to update the solution at the boundaries. At r = L, our 
conditions (2.24), (2.26), (2.27) are all of the form: 

^ = <?(/>,*)• (47) 

A second order discretization of this is given by: 

S(U£; At ) = S(tfr +1 ) + S(U' N ) - 5(C# A£ ) + 2AtQ(U^). (4.8) 

Note that all quantities on the right are available from the interior scheme. 
Another numerical condition is needed. We obtain it from the equation for 
the outgoing characteristic, (2.4). Writing it in the form: 

dR . dR . . 

we have the second order discretization: 

(1 + C) R{V$g) = *(ffjt +1 ) + R{U l N ) - R(Uj^ At ) + W 

(4.9) 

- C(R(Uh +i ) - R(U' n ) - R(U' n +a % 


where 



and 

W = 2A tW(u£ + J). 

These equations yield updates of the Riemann variables at the artificial 
boundary. Equations (2.1) and (2.2) are inverted to update the primitive 
variables: 

p(r N+1) t + At) = G _1 (^ 2 ^) ’ (4.10) 

z{r N+1 ,t + At) = ^p{r N+1 ,t + At){R + S). (4.11) 

At the origin, continuity of the velocity field requires: 

z(0,t + At) = 0. (4.12) 


Again, a numerical boundary condition is required which we obtain from 
the characteristic equation, (2.5). Writing it in the form: 


as N ds 


W(p,z), 


we use the second order approximation, 

(l - d) S{U [ +At ) = S{U{) + S{U\) - 5(tf 2 t+A< ) + W 

- C(S(U*) - S{U\) + S(tf 2 t+At )), 


where 


and 


5(1/1+**) At 
2 

Ar 


n 1, 


W = 2AtW{U l 3 ). 

2 

Given 5, p can be computed: 

p(0,t 4- At) = G-\-S ). 


(4.13) 


(4.14) 


This completes the update of U i 


t+&t 
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NUMERICAL EXPERIMENTS 


Here we present some numerical calculations to validate the effectiveness 
of our boundary conditions. For purposes of comparison we include the 
condition resulting from Thompson [10]. We here list and label the different 
conditions that we have used for the computation. All of them have the 
form 

I G - cw ) - c (51) 


where 
(HH1) Q 

(HH2) Q 

(HH3) Q 

(Th) Q 


*%/ /'(poo) 

pL 

~ <?(*»)) 

(*(£,*) - G( Poo )) 

, / # (Poo ) 

L 


(from equation(2.26)) 
(from equation(2.27)) 
(from equation(2.24)) 
(from Thomps on [ 7] ) 


A simple idealized explosion problem is considered in which the density 


is initialized to 

Po{r) = 

and the momentum is initialized to 


3 r < 1 
1 r > 1 


*o(r) 


0 0 < r < L . 


The true solution includes the propagation of a weak, decaying shock 
with the solution on the truncated region eventually approaching the steady 
state p = 1 and z = 0. We note that only this steady state is compatible 
with boundary conditions HH2 and HH3. By contrast, HH1 and Th are 
compatible with steady states at any density. In the graphs presented here 



the momentum and density are plotted against r at different time steps. We 
use a spatial mesh width of A r = .05 and a time step chosen so that the CFL 
number corresponding to sound speed of the compressed gas, yj f f \ 3), is .25. 
In the first case considered here the far field boundary is located at 5 (L = 5). 
Figures 1-4 (a) show, at time steps 400, 600, 800 and 1000 respectively, the 
results obtained from the computations using the condition HHl. Similarly 
for conditions HH2, HH3 and Th the results are reported in figures 1-4 (b)- 
(d). As can be seen in the figures the solutions are initially qualitatively 
the same for each boundary condition. For longer times, however, marked 
differences in the solutions appear. All approached a steady state. As 
discussed above, this is necessarily the correct steady state for HH2 and 
HH3. For HHl the final density was roughly .993, an error of about .7 %. 
For Th it was .984, an error of 1.6 %. 

The contrasting results are accentuated by further contraction of the 
computational domain. Figures 5 and 6 show the results obtained for L = 
2.5 employing boundary conditions HH3 and Th respectively. Even at time 
400, the results obtained with the nonreflecting condition are seen to be 
significantly in error, while those obtained with our asymptotic condition 
are not. Again, the steady state density found using HH3 is correct while 
that found using Th is off by about 12.5 %. 

It is worthwhile to note the significant fluctuations in the variables which 
occur near the origin. We believe this is a natural consequence of the fo- 
cussing of incoming spherical waves. The diminished amplitude of these 
fluctuations resulting from Th on the smallest domain is evidence of its in- 
accurate representation of the transient solution. That is, some physical, 
incoming waves were not generated. 

Computations were performed on a Sun Microsystem 3/260 with a float- 



ing point accelerator. It took about 56 seconds of cpu time for the first case 
( L — 5) and 30 seconds for the second case ( L = 2.5). In each case the total 
number of time steps was 2000. 

In conclusion, we have established the accuracy of our reflecting condi- 
tion even when implemented on a domain of modest size. Such ideas become 
crucial in truly multidimensional computations. Though we have not per- 
formed any such computations, we present ideas in the next section for the 
generalization of our boundary conditions to nonsymmetric problems. 

GENERALIZATIONS TO NONSYMMETRIC FLOWS 


We now consider the Euler equations in spherical coordinates for non- 
syminetric, isentropic flows. The new dependent variables are taken to be 
the angular momenta, m and q : 


dp dz 1 dm 1 dq __ 

dt ^ dr r d6 r sin 6 d<j) 


(6.1) 


dz d ( z 2 . \ 13 ( mz\ 

at + Tr\J + M ) + rTeWJ 


1 d fqz'' 


dm d ( mz\\ 


i d_ 

' + rde 


m 


+ /(/>) + 


dt dr \ p J 

d fqz\ 1 d / qm\ 1 

dr \ p ) r dO V p ) r sii 


r sin 9 d(j> \ P 
1 d ( mq\ 

r sin 9 d<j> V p J 


dq 

dt 


sin 6 d<j> 


\ + fiP) 

t P i 


where we have introduced, 


= 92 , (6-2) 
= <73, (6.3) 
= <74, (6.4) 


9 1 


2 z m cot 9 

_ , 

r r 


92 = 
93 


2 z 2 — m 2 — q 2 mz cot 6 


pr 


pr 


3 mz cot 6. 9 9. 

-(m 2 -q 2 ) y 


pr 


pr 
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9 a = 


3 qz 
pr 


2<?racot 9 
pr 


We again consider a domain exterior to a compact body and assume the 
initial conditions satisfy p — p^ and z = ra = g = Oforr>T. Following the 
construction of section 2, we work with the symmetric Riemann variables, R 
and S', as well as the angular momenta. The formal expansion we postulate 
has the form: 

*(r, >,M = Ho + MjflM 


S(r, 0,^,1) = 5o + 5l ^’ - — + 

V 


m(r ,$,</>, t) = 
q{r,0,4> t t) = 


r 2 

qA r ilA 1 t l 


+ 


r 2 + 

(6.5) 

S 2 {r,e,<f>,t) 

7 * ' * • } 

T l 

(6.6) 

m 3 {r, 0, <j>, t ) 

o *T • * * > 

T 6 

(6.7) 

qa{r,6,<l>>t) 

n 1 ’ * * ? 

r J 

(6.8) 


where, borrowing results from the linear case, we assume the angular mo- 
menta, m and g, are O(^) as r — > oo. The equations for i?i, i? 2 ) $i and 
^2 are taken unchanged from section 2. This involves the neglect of lower 
order terms involving 6 and <f> derivatives. Although not entirely consistent 
with our inclusion of lower order terms involving ^ , this is justified by the 
expectation that the primary direction of propagation in the far field is the 
radial one. Expressions for the first corrections may then be copied from 
(2.14), (2.16) and (2.17): 


<(r,0,^,r) = r + 


Si{r,0,<t>,t) = 0, 

Ri{r,O,<f>,t{r,0,<f>,r)) = 

6 t) In [ %$»$ ] 
s/TTp oo) 


(6.9) 

( 6 . 10 ) 

( 6 . 11 ) 
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More importantly, the boundary conditions (2.24), (2.26) and (2.27) are 
unchanged. Equations for m 2 and q 2 are: 


Here, 


dm 2 

dt 


21 dm 2 dpi n 

+ 5 — + / {Poo)-^r = 0 , 

r Poo dr 60 


d$2 Jl_0q2 f{Pno)dpi 
dt rpoo dr sin# d<j) 


z i 


PooRl PooH 

2 ~ 2 


( 6 . 12 ) 

(6.13) 


and 


n ^00^1 Poo B 

pl ~ 2 7Kp~) ~ 2y/r( Poo y 

We include the * term so that characteristics can be computed. In partic- 
ular, at a point of outflow, z x > 0, no boundary condition for m 2 and q 2 
is required. At inflow we may simply use (6.12) and (6.13) without the r 
derivative terms to update the angular momenta. In summary we have: 




2L 

Lp(L 1 9 1 4> 1 t) 

L 


and, if z(L 1 6 y <f) i t) < 0, 




(6.14) 


(6.15) 

(6.16) 
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16. Abstract 

The numerical solution of exterior problems Is typically accomplished by 
introducing an artificial, far field boundary and solving the equations on a 
truncated domain. For hyperbolic systems, boundary conditions at this boundary 
are often derived by imposing a principle of no reflection. However, waves with 
spherical symmetry in gas dynamics satisfy equations where Incoming and outgoing 
Riemann variables are coupled. This suggests that 'natural 1 reflections may be 
Important. We propose a reflecting boundary condition based on an asymptotic 
solution of the far field equations. We obtain nonlinear energy estimates for 
the truncated problem and present numerical experiments to validate our theory. 
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